/*
 * StructureFactor.cpp
 *
 *  Created on: 29.07.2012
 *      Author: maikel
 */

#include "StructureFactor.h"

#include <fstream>

#include "Atom.h"
#include "Timeset.h"
#include "Eigen/Dense"

using namespace std;

namespace gro {

void StructureFactor::skip_line_until_valid_tilte() {
	while (!timestream.eof()) {
		string line;
		streampos start = timestream.tellg();
		getline(timestream, line);
		if (line.find("t=") != string::npos) {
			timestream.clear();
			timestream.seekg(start);
			break;
		}
	}
}

StructureFactor::StructureFactor(string grofile) : N(0) {
//	timestream.exceptions(ifstream::failbit | ifstream::badbit);
//	recovery.exceptions(ofstream::failbit | ofstream::badbit);
	timestream.open(grofile.c_str());
	recovery.open("last_calc", ofstream::trunc);
}


double StructureFactor::sum(Vector3d q, const vector<Atom> &atoms) {
	double sum = 0;
	this->N = atoms.size();
	for (size_t i = 0; i < atoms.size(); i++) {
		for (size_t j = 0; j < i; j++) {
			sum += cos(q.dot(atoms[i].position - atoms[j].position));
		}
	}
	return sum;
}

double avg(vector<double> x) {
	double sum = 0;
	size_t i;
	for (i = 0; i < x.size(); i++) {
		sum += x[i];
	}
	return sum / x.size();
}

double StructureFactor::operator()(Vector3d q) {
	if (cache.find(q) != cache.end()) {
		return cache[q];
	}
	streampos start = timestream.tellg();
	vector<double> sums;
	Timeset timeset;
	while (!timestream.eof()) {
		try {
			timestream >> timeset;
			double s = sum(q, timeset.atoms);
			sums.push_back(s);
			recovery << timeset.time << " " << s << endl;
		} catch (Timeset::NoTime &error) {
			skip_line_until_valid_tilte();
		}
	}
	return 1 + 2 * avg(sums) / N;
}

} /* namespace gro */
